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FOREWORD 


This report covers the third year of a research effort 
devoted to the determination of minimum noise aircraft land- 
ing trajectories. Increased concern for environmental pro- 
tection, as well as improved measurement and instrumentation 
capabilities, have provided the primary impetus for this 
work. Our study has been concerned with the Boeing 737, a 
short-haul passenger aircraft, and the Patrick Henry Airport 
which is located at Newport News. 

During the three years for which this research has been 
in progress, we have employed in addition to the principal 
investigator one post-doctoral researcher and one masters 
candidate who just received his degree. Also, another mas- 
ters student has been employed part-time during one summer. 

Besides the three annual reports, there have been two 
technical papers written on our work. One of these has been 
published and the second one is being reviewed. 

Our goal is a working computerized optimization program 
which may be modified by changing the population data to 
yield optimal trajectories for any airport. This report 
gives the current status of the effort. 
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I. INTRODUCTION 

The objective of this report is to bring the reader up 
to date on the present status of our landing trajectory opti- 
mization research without repeating a great deal of what has 
already been reported in the reports [1,2] on the effort 
during the first two years. 

The accomplishments during the first two years included 
the following: 

1. Develop the aircraft equations. 

2. Adapt wind-tunnel data for computer usage. 

3. Obtain a passenger comfort model. 

4. Develop a noise model. 

5. Develop a population model. 

6 . Integrate the noise model and population model . 

7. Establish a perf03onance measure. 

8. Use the performance measure to compare constant 
glide slope trajectories. 

9. Set up the equations for the steepest descent 
optimization procedure. 

Each of these items, except number 2, has been discussed 
in detail in the two previous annual reports. For complete- 
ness, discussion of this item will be included here. In 
addition, the report covers the accomplishments of the past 
year which are the following: 

1. Programming and modifying the steepest descent 
optimization procedure. 
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Successfully iterating toward the optimum for a 
four-mile trajectory. 

3. Beginning optimization runs for a twenty-mile 
trajectory. 

For reference the two technical papers which we have 
written are included in the bibliography [3,4]. 



II. GENERAL SYSTEM EQUATIONS 


2 . 1 Introduction 

For the discussion on our treatment of the wind tunnel 
data to be meaningful, one must have in mind how these data 
are to be used. Thus, a brief description of our simulation, 
along with definitions of the variables, seems in order. The 
realistic simulation of aircraft behavior necessitates the 
solution to the nonlinear, differential equations of motion. 
These equations have been formulated as first-order deriva- 
tives of the state variables, and they describe the complete 
six degrees-of-freedom of an actual aircraft. 

The equations of motion can have many forms; the speci- 
fic form depends upon the choice of the coordinate system. 
Figure 2.1 illustrates the relative orientation between three 
■possible systems. The origin of each of the systems is the 
aircraft center of gravity. 

The body axes (X,Y,Z) are rigidly fixed to the aircraft. 
The Y axis is perpendicular to the aircraft's plane of symme- 
try and is directed out the right wing. The X axis is in the 
plane of symmetry and points toward the front of the air- 
craft. The Z axis is normal to the X-Y plane and forms a 
right-handed system. The fuselage reference line (FRL) coin- 
cides with the X axis. 

The stability axes differ from the body axes by the 
angle of attack (a_„_). The X axis lies in the plane 

J! XUj S 

defined by the relative wind vector and the Y body axis, the 

latter coinciding with the Y axis. The Z axis is 

s s 
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X,Y,Z - Body Axes 

X ,Z - Wind Axes 
w vr w 

X^,Y ,Z - Stability Axes 
s s s 

Axes Orig.in at the Center 
of Gravity 


Figiu'C- 2,1 


Coordinate System Relationships 
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perpendicular to the X -Y plane and forms a right-handed 

o S 

system. 

The wind axes have the X axis collinear with the rela- 

w 

tive wind vector and differs from the X axis by the side- 

s 

slip angle (p). The axis co^-ncides with the axis and 

Y is normal to the X -Z plane, 
w w w 

The body axes representation yields simple expressions 
for the Euler angles of the aircraft. -However, the relative 
wind vector, glide slope, angle of attack, and side-slip 
angle are more difficult to calculate. The wind axes permit 
simple calculation of the translational equations, angle of 
attack, and side-slip angle. The disadvantage lies in the 
complexity of the moment equations and the variable inertia 
values . 

To minimize the complexity of the model, a combination 
of wind and body axes is used [5]. The translational equa- 
tions are solved in the wind axes, and the rotational equa- 
tions are solved in the body axes. The flight path coordi- 
nates require the transformation from the body to the earth- 
fixed axes. 

The earth-fixed axes have the X^, and Y^ axes mutually 
perpendicular, located in the ground plane. The Z_ axis is 
directed normally into the ground. The center of this system 
is dependent upon the population model# 

The aerodynamic data for the aircraft being modeled was 
obtained using the stability axes system. This system per- 
mits the simple derivation of the aerodynamic forces and 



Figure 2,2 
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moments. These forces and moments must be resolved into the 
appropriate system for use in the equations of motion. 

Figure 2.2 is a flowchart of the sequence used to solve 
these equations. This sequence is outlined in the follow- 
ing section. 

2.2) System Equations [6,7^ 


The sequence begins by expressing the thrust vector 
from each engine in terms of stability coordinates. 

T = F + F 


T 

s 

n n 

= i 

X 

T 

=■ 0 

s 

y 

T 

= -T sinO-Ti^T-, 

s 

z 

WCP 

“■WCP 

= °^FRL **■ 1 


(2ol) 


The thrust vectors were assumed to be parallel and of equal 
magnitude. The angle of attach of the engines is the same 
as that of the wing chord plane (WCP) , The moments induced 
by the thrust vectors are determined in the body axes. 

I+P = = 0 

= T (2.2) 

The load factor is also evaluated, 

^ (2.3) 

z 

Next, the gravity forces are resolved into components 
in the stability axes. 
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Fg = W ( cos0cos0sinap^^ - sinecosapj^^^) 
X 

Fg = W cos9sin0 

y 

Fg = W ( sinesinttpj^^ + cos0cos0cosapj^j^) 
z 


(2.4) 


The angular velocity compqnents of the aircraft are 
then expressed in the stability axes. 


Pg ~ P ■** ^ ^^^“•FRL 


qg = q 


(2.5) 


r^ = -p sinapj^L + 


cosa 


FRL 


The aerodynamic forces and moments are calculated in 
the stability axes system. These forces are then combined 
■with the thrust and gravity forces and transformed into the 
■wind axes system. 

F =r ( T + Fp + X ) cos3 + ( Fp + Y ) siriP 

W S S Lr S 

X XX y 

y y XX 


( 2 . 6 ) 


F = T + F^ + Z 

s 

z z z 


This is followed by the integration of the velocity, angle 
of attach, and sideslip angle. 


= V = F /m 


f^ = a = ( -qy'cosp) - pg tanP + 


(2.7) 


■^';here 


^3 = ^ 


q = -F /mV 
w w ' 

z 

r = F /mV 
w Wy' 
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The aerodynamic moments are then transformed to' the body 
axes system. 


L = cosapj^j^ - Ng sinapj^j^ 

^ ^ ^s (2.8) 

N = sinapj^^ - cosa^^^ 

The usual form of the roll, pitch, and yaw moment 
equations is - 

^xx ^4 ~ ^xx P “ ^^yy ■ ^XZ 

+ L + 

lyy = lyy q = - i^) rp + (r^ - p 2 ) 

+ M + Hp (2.9) 

^zz ^6 “ ^zz ^ ” ^^xx “ ^yy^ ■*■ ^xz - *5^^ 

+ N ^ N^. 

Equation (2.9a) contains an undesirable derivative (r) in 
the right hand side. To eliminate this derivative, equa- 
tion (c) was substituted into (a). The final form of the 
. equations is 


^^xx^zz " ^iz^ ^4 “ ^^yy^zz “ ^Jz " 

**■ ^^zz ■*' ^‘xx-“ ^yy^ ^xz^^^ 

^yy ^5 "" ^^zz ~ ^xx^ ■*' -^xz^^^ - P^) + M + 

( 2 . 10 ) 

^ZZ ^6 ~ ^^xx “ ^yy^ P'3 ^xz ^P “ '3^) + N + . 

Next, the Euler angular rates are integrated. 

f 7 = V = (r COS0 + q sin0)/cos9 • 

» 

fg = 9 = q COS0 - r sin0 


( 2 . 11 ) 
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f .= = p + f sine 

9 7 

Integration of the aircraft velocity components in the 
earth-fixed axes is the next step. 

f^Q = Xg = V(cosapj^j^cos3cosecosy^ + sinp (-cos^siny' 
+sin0sin0cosj(;)+ sinap^^cosp (sin0sinjy 
+ cos0sin9cos(//) ) 

= Yg = V(cosapjy^cospcos 0 sin^ + sin? (cos 0 cos 51 / 

+ sin 0 sin 0 sin^) + sinapj^^cos? (~sin 0 cos'p 
+ cos 0 sin 0 sin^) ) . ( 2 . 12 ) 

f ^2 “ “ y(cosappj^cos?sin 0 - sin? (sin0cos9 ) 

- sinapj^j^cos?cos9cos0 ) ) 

An additional differential equation vas included into 

the system for the performance index. It has the form 

« • 

= J = K^Wp+ KgP K 3 + PENALTY FUNCTIONS. (2.13) 
where the K^’s are scaling constants, is the fuel con- 
sumption rate, and P is the instantaneous population 
•exposed to noise. The constant is the constant for the 
time component of the performance index. 

2.3) Integration Method 

2 . 3 . 1 ) Runge-Kutta Fourth Order 

An examination of the differential equations (2.7), 
(2,10), (2.11), (2.12), and (2,13) reveals that they are 
highly nonlinear and require numerical integration. The 
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original method employed was a Runge-Kutta fourth order 
algorithm. This algorithm uses a weighted average of four 
estimates of the dependent variable over the interval At, 
to obtain the value at t+at. 
ki = f(x^^,t^)At 
Rs = f (x^+Xi/2,t^+At/2)At 

^3 = f (x^+k2/2,t^+At/2)At (2.14) 

)C4 = f (x^+k3,t^+At )At 

x(t+At) = (Ki+2Kl3+2k3+R4 )/6 + x(t) 

The reliability and self -starting characteristics of this 
method were the motivating factors in the original selec- 
tion of this algorithm. The choice of a suitable integ- 
ration interval still remained, A trial and error procedure 
was employed to determine the appropriate step size. A 

of simulations were examined, each with a success- 
ively larger step interval. When two consecutive simulations 
differed appreciably in their state histories, the next 
®tsp interval tested was the average of the previous two 
values. This procedure resulted in the final choice of 0.1 • 

seconds as the step interval. The ideal choice would have 
been an interval which satisfied the relation A.t= 2“^, where 
n is an integer. Excessive computer execution time preven- 
ted the use of At=0,0625 seconds (n=4); and At=0.125 sec- 
onds (n=3) did not yield sufficient accuracy. With an 
interval of .10 seconds, the Runge-Kutta algorithm still 
consumed an undesirable amount of execution time. This 
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supplied the motivation to examine more rapid integration 
algorithms. 

2.3.2) Milne-Reynold’s Method 


The method finally chosen was the Milne-Reynold’s 
second order predictor -corrector algorithm [[7 ] . This 
technique requires values of the dependent variables .. . • 
at the four previous discrete time points and is not self- 
starting. The Runge-Kutta algorithm is utilized for the 
first three time intervals to supply the three additional 
points. The Milne-Reynold’s method is then used to com- 
plete the simulation. 

The algorithm first estimates the derivative using the 
current (t=t^) state variables. A predicted value of the 
state at t=t^+/^t is then obtained from the relation 


-f ,+ 2 f 

n+1 n-3 3 n n-1 n-2 


(2.15) 


These predicted values, used in the system equations (2.7), 
( 2 . 10 ), ( 2 . 11 ), (2.12), and (2.13), provide a revised esti- 
mate of the state variable derivatives. 

(2.16) 

The corrected values for the state variables at t +^t are 

o 

then calculated . 




(2.17) 
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The 0.10 second step interval was retained and a comparison 
of tliG Runge— Kutta and Milne— Reynold ' s methods was performed. 
Variations from the Runge-Kutta algorithm of 10~"^ percent 
in displacement and 0,10 percent in the pitch rate were 
observed for a 300 second trajectory. Such a deviation is 
considered acceptable. 

The reduction in the computer execution time was con- 
For each second of sxmulated txmej the Runge- 
Kutta method required 1.07 seconds of execution time, the 
Milne-Reynold’s method only 0.61 seconds. Thus, a 43 
Percent reductxon xn executxon txme was achieved by using 
the latter method. 

2 . 4) Summary 

The equations describing the aircraft.' s motion have 
been programmed. Initially a fourth order Runge-Kutta 
algorithm was considered as the integration method. , However, 
the execution time proved excessive and a more suitable 
algorxthm was sought. The algorithm finally chosen was the 
Milne-Reynold's predictor-corrector method. This method 
was chosen for its speed and stability characteristics. As 
a second order technique it required 43 percent less execu- 
tion than did the Runge-Kutta. For every simulated second 
of flight the Milne-Reynold's method required 0.61 execution 


seconds . 



3,0) Aircraft Model 


3.1) Introduction 

The nonlinear differential equations given in Chapter 
2 are the. general relationships governing the motion of an 
aircraft. The specific aircraft characteristics for a 
Boeing 737 are manifested in the dependence of the aerody- 
namic forces (Xg,Yg,Zg) and moments (Lg,Mg,Ng) on the 
vehicle state and control variables. . An accurate model of 
these characteristics is necessary if the optimal trajec- 
tories are to be physically realized. In addition to the 
~ modeling of the aerodynamic data, the physical constraints 
on the aircraft must be faithfully reproduced. Constraints 
such as descent and ascent rates, maximum altitude, and 
speed restrict the set of trajectories from which the 
optimum can be determined. The following sections of this 
chapter detail the evolution of an appropriate aircraft 
model for a Boeing 737-100 utilizing two JT8D-7 turbofan 
engines . 

3.2) Source of Data 

Figure 3,1 illustrates the configuration of a Boeing 

737-1-00 aircraft and includes the pertinent dimensions [S"]. 

Aerodynamic data for this aircraft was obtained in the form 

of stability derivatives, supplied by the Flight’ Instrumen- 

ts 
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tation Division personnel at NASA*s Langley Research Center* 
Figure 3*2 is a typical graph and illustrates the .dependence 
of the basic- lift coefficient on the angle of attack. and the 
setting. The equations which define the total aerody- 
namic coefficients are presented in the Appendix, 

3,3) Stability Derivative Simulation 

The major task involved in developing the aircraft model 
was to adapt the aerodynamic data into a form more amenable 
to digital computation. The two alternate methods of adap- 
tation considered were the direct storage of the graphical 
data and a functional approximation. Each has advantages 
and disadvantages which are described in the following sec- 
tions . 

3,3,1) Direct Data Storage 


This method is a piecewise linear approximation of the 
data. For example, in Figure 3,2 each curve would be 
represented by data points chosen from that curve. The 
choice of the specific points depends upon the degree of 
nonlinearity, A second influential factor is the choice of 
either uniform or nonuniform abscissa intervals. The esti- 
mated computer storage requirement for uniform intervals is 
3,000 decimal words. The abscissa intervals and ranges 
assumed in this estimate are given in Table I, 
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Figure 3,2 


Basic Lift Coefficient 
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Table I 

Abscissa Uniform Interval Sizes and Ranges Used in 
the Direct Data Computer Storage Estimate 


Abscissa 

Variable 

Step 

Size 

Total 

Range 

'Angle of Attack 

3.0° 

-5°' to 5° 

Mach Number 

0.2 

0 to 0 , 6 

Rudder Deflection 

5.0° 

-25° to 25° 

Aileron Deflection 

10.0° 

0° to 40° 

Flight Spoiler 
Deflection 

0 

o 

• 

in 

0° to 40° 

Sideslip Angle 

o 

o 

• 

in 

-15° to 15° 


The use of nonuniform intervals required an estimated 
total of 1,600 data points. This estimate required fewer 
points per curve than the uniform interval method to achieve 
the same accuracy. However, a knowledge of both the ordinate 
and abscissa at each point was mandatory. The computer 
storage requirement was 3,200 decimal words. The logic 
needed for implementing either method was nearly equal. 
Therefore, if direct data storage were chosen as the adapta- 
tion technique, then uniform abscissa intervals would be 
more efficient, 

3.2.2) Functional Approximation 


The. second alternative involved the identification of 
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functions which approximate the data to an acceptable degree 
of accuracy. A least-squares criterion was. employed to 
generate the desired functions, A polynomial representa- 
tion was chosen as the nominal function. 

Figure 3.3 illustrates the comparison of two curves 
from Figure 3.2 and their respective polynomial approxima- 
tions. These approximations are typical of the accuracy 
obtained and were judged as acceptable. The equations for 
the polynomials in Figure 3.3 are 

= 1,218 X 10“^ t- .1072 (3.1a) 

-3 ■ -4 

+ 3.266 X 10 " ^‘247 X 10 

+ 1.173 X 10“^ 


for a flap setting of zero and 


Base 


1.216 + ,1742a 

6.353 X 10”^ a 


WCP 

WCP 


3.306 X 10 ^ a 


WCP- 

(3.1b) 


for a flap setting of 40, 

The computer program, which implemented the least - 
squares algorithm, fitted successively higher order poly- 
nomials to the data. The program terminated when the 
polynomial estimated the data points within an error of 
two percent . 

Whenever possible, the coefficients of similar order 
polynomials were themselves fitted with least-squares . ■ 
polynomials in an effort to reduce the computer storage 
requirements. As an example. Figure 3.4 shows curves for 
each of five values of altitude (0, 13000, 20000, 23000, 



^’keceding page 



Mach Number 


Figure 3,4 

Variation of the Pitching Moment Coefficient 
For Changes in the Pitch Rate 


N.' 



and 35000 feet) . 
polynomials 


These curves were approximated by the 


(— 

dq 


.25 


) 


/ = 0.2533 - 0.02075 M 

I = 0.2519 - 0.00286 M 

= 0.2522 - 0.00554 M 

= 0.2518 t 0.00839 M 

^ = 0.2508 + 0.01870 M 


These equations have the form 


Y = C + C, M. 
o 1 


h = 0 feet 
h = 13000 
h = 20000 (3.2) 

h = '23000 
h = 35000. 


(3.3) 


By. treating each as a function of the altitude, a least- 
squares approximation for each coefficient was obtained. 


The result was that the five curves of Figure 3.4, which 
originally required the ten coefficients of equations (3.2), 
were represented by the single equation 


dC. 


M 


.25. 


dq 


(D^ + D^h) + (D 2 + D^h + D^h2) M; 

h = (0, 13000, 20000, 23000, 


(3.4) 

35000 feet) 


where • 

D = 0.2532 
o 

•D^ = -6.671 X 10“® 

= -2.078 X 10 “^ 

= 1.542 X 10"^ 

= -1.177 X 10"^^, 

and which contained only five coefficients. 

For many stability derivatives the two percent .error 
criterion required a polynomial of an undesirably high order. 
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For these derivatives piecewise polynomials of fourth order 
or less were used. Figure 3,5 shows a typical curve for 
which this method was applied. The polynomials and their 
applicable regions are 


d^ 


0.6582 + 8.796 x lO'^ 

2.412 X 10“^ agcp + 2.222 X 10“^ 

2.689 X 10-® ^ <^WCP ^ 20°) 

1.826 + 0.1458 a^cP " 2.973 x 10“2 

(20° < 25°). (3.5) 


In a further effort to reduce the number of coefficients, 
Fourier cosine series were fitted to the piecewise polynomials. 
These series approximations were originally used in twelve 
instances. The form of the series is 

Y = Cq + I COs[(-^p — ) + (3.6) 

k=j. 

If the series required fewer coefficients (C^ and 0 ^) then 
the latter representation was used, otherwise the polynomials 
•were retained. A second motivation to employ a Fourier 
series was the presence of roundoff errors. If the inde- 
pendent variable in a polynomial was the aircraft altitude, 
the difference between two nearly equal, large numbers was 
obscured by the roundoff errors. ' 

Subsequent examination of the twelve Fourier series 
indicated that the accuracy of the approximation was not 
sufficient for seven of the derivatives. The piecewise- 




Angle of Attack (deg) 


Figure 3«5 

Piecewise Fitted Polynomials 
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polynomials were reinstated for these stability derivatives. 
The functional approximation technique required 2,000 
decimal words of computer storage. 

3.3.3) Choice of Representation 

The software logic complexity required for either direct 
data storage or functional approximation was estimated to be 
equal. Based upon the difference in the storage require- 
ments of the two alternative methods, the functional approx- 
imation model was chosen. The effort to model all the sta- 
bility-derivatives by polynomials required an estimated 
200 man-hours. An additional 200 man-hours were needed to 
develop an error-free software model. If computer storage 
is not a primary concern, considerable time could be saved 
by the use of direct data storage. The relatively short 
time to implement this data representation (20 to 40 Man- 
hours) greatly facilitates modification of the entire soft- 
ware package for an aircraft other than a Boeing 737-100. 

3. -3. 4) Implementation 

Returning to Figure 3.2, there are nine separate 
curves which describe the variation of the basic lift coef- 
ficient as the, angle of attack changes. The choice of the 
appropriate curve is governed by the current setting of the 
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pilot flap quadrant located in the cockpit. This quadrant 
has nine disrete settings (0,1,2,5,10,15,25,30,40). Consi- - 
derable programming complexity would have been required if 
the simulation permitted only these discrete flap settings, 

A decision was made to permit a continuous variation in flap 
setting in the range 0 to 40, The evaluation' of any stability' 
derivative at the intermediate flap settings is achieved by 
linear interpolation. 

To illustrate the technique, assume the basic lift 
coefficient of Figure 3,2 is desired for an angle of attack 
of 5°' and a flap setting of 32, First, the coefficient 
values for_ apjQp=5° and flaps at 30 and 40 are calculated to 
be lo59 and 1,96, respectively. The value for a flap set- 
ting of 32 would then be interpolated to be 1.66, This 
same technique was anployed for all curves having multiple 
parameter values , The variables used as parameters in modeling 
the aerodynamic data were flaps, rudder deflection, altitude, 
and angle of attack. 

Special attention was directed to the modeling of the 
effectiveness factors for the flight spoilers and ailerons, 
.The flight spoiler panels 2&3 can be deployed inde- 

pendent of panels 6&7 ^^ 5 gp 2 ^« Therefore, separate values 
of ^6sp calculated for each set of panels. Reference 

to Figures 3,1 and 3,6 indicate that both sets of spoilers 
produce negative lift and positive drag and pitch. Their 6on- 
tributions to the side force, roll and yaw moments reguire • 
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'^6spl"=° ailerons are deployed symmetri- 
cally. A positive value for (aileron effectiveness 

factor) occurs when the right aileron trailing edge is up. 

3.4) Model Restrictions 

A number of additional assumptions were made in the 
development of the aircraft model. Since the trajectories 
of interest are those of the landing phase only, the Mach 
number is not expected to exceed 0.6. Thus all stability 
derivatives dependent on the Mach number have been modeled 
only for speeds below a Mach number of 0.6. All trajec- 
tories should be examined to insure that the Mach number 
does not violate this assumption. If such a condition 
does occur, significant deviations from actual vehicle ' 
performance will result. 

The ground effects in the aerodynamic data have been 
ignored. Thus the minimum "relative"^ altitude for which 
the aircraft model is valid is 100 feet. Below this alti- 
tude ground effects become significant. 

The ground spoilers on the aircraft were assumed to 
have been locked in the undeployed position. These controls 
are used' only as speed brakes, when the aircraft is decelera- 
ting during the ground run phase of the landing maneuver. 

^"relative" altitude is defined as the height of the 
aircraft above the local ground plane. 
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The restriction on the minimum altitude eliminated- the need 
to model these control surfaces. 

Throughout the .program the altitude is evaluated rela- 
tive to the mean sea level. Thus the initial aircraft 
altitude for a trajectory is the sum of the airport altitude 
and the desired height of the vehicle above the local 
ground plane. The mean sea level altitude of any airport 
can be obtained from a document similar to £ 9 3* 

The atmosphere used in the simulation is from [lO]. 

The variation of the air density and speed of sound with 
altitude are described by (3.7). 

(e = 2.38 X 10"^ (1 - 0.668 

-5 (3.7) 

a = 1117 exp( -0.36 x 10 h) 

In addition, the atmosphere was assumed to have no wind 

velocity. 

The aircraft total weight and moments of inertia were 
•assumed constant over the entire trajectory. The effect of' 
the constant weight assumption was tested in simulations 
restricted to the longitudinal plane. The variations in 
.the aircraft performance are summarized in Table II. The 
largest deviation for a 180 second simulation does not 
exceed three percent. The moments of inertia were assumed 
constant for the tests. The variations were judged to be - 
acceptable . 

The fuel consumption of the JT8D-7 engine is, nominally. 
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Table II 


Effect on Performance of The Constant 
Weight and Moment of Inertia 
Assumption 


State Variable Variable Constant Variation 

3t t^= 180 (s) Weight Weight {%) 


Angle of Attack: 
tdeg) 

3.7896 • 

- 3.7801 

0.25 

Velocity 

(ft/sec) 

209.02305 

209.39176 

-0.18 

Pitch Rate 
(rad/sec) 

-8.353 X 10"^ 

-8.4375 X 10"^ 

-1.01 

Pitch Angle 
(deg) 

1.01496 

0.98450 

3.00 

iX Groimd 
Coordinate 
(ft) 

-7206.167 

-7199.2207 

C.IO 

Y Ground 
Coordinate 
(ft) 

29537.837 

29505.158 

0.11 

Altitude 

(ft) 

1686.0828 

1675.8523 

0.61 
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0.. 6 Ibs/hr/lbf , For a thrust of 7,000 pounds, this is a 
fuel flowrate of 2.33 Ibs/sec. If the total aircraft is 
90,000 pounds and the trajectory duration is 500 seconds, 
the fuel consumed is 1,3 percent of the original weight. 

The added complexity of simulating the inertial time depen- 
dence was not considered justified for such a small increase 
in accuracy. 

The landing gear was assiamed locked in the deployed 
position during the entire trajectory. This assumption is 
not considered critical, since the contribution of the 
landing gear to the aerodynamic forces and moments is rela- 
tively small. In addition, the elevator was considered inde- 
pendent of the stabilizer. This was not strictly true and 
can be modified by specification of a control constraint. 

A few of the stability derivatives were modeled in a 

special manner. For derivatives such as (C ),y(C ) 

nsp'^M'^ ' nsp^M=0’ 

which have as the parameter, the value for is 

evaluated assuming a^^p=0. Similarly, when the 

stability derivative is calculated asstiming Another 

derivative which required special consideration was (dC^/dp). 
This derivative is given for both high (M > 0.4) and low 
speeds. In an actual flight, the flaps are not deployed 
for high speeds. Under this condition (dC^^/dP) has the Mach 
number as the parameter. For low speeds the parameter is 
flap setting.. Due to the open loop nature of the control 
vector, it is possible for the flaps to be deployed when the 
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Mach number exceeds 0.4. When this occurs, the program 
evaluates the stability derivative as if the Mach number 
is less than 0.4, The use of penalty functions in the 

i 

performance index can aid in removing these unrealistic 
conditions. Specifically, penalty functions which enforce 
the stall, buffet, and load factor constraints will heavily 
penalize the performance index whenever the flaps are 
deployed for Mach numbers in excess of 0.4. The optimi- 
zation algorithm would then minmize the performance index 
by specifying a control history which would avoid these 
unrealistic conditions, 

3.5) Aircraft State Constraints 

If the simulation is to faithfully model the aircraft ,- 
the physical constraints on the vehicle must be enforced. 
This was accomplished by the addition of penalty functions 
to the performance index. The various constraints and 
■^heir respective penalty functions, which must be satisfied, 
are described in the paragraphs . which follow. 

3.5l) Maximum Ascent and Descent Rates 

The maximum rated ascent (descent) rate for the 737- is 
100 (250) ft/sec. Whenever the glide slope is positive, 
the aircraft is climbing and the ascent rate is compared to 
the specified maximum. Therefore, the penalty functions have 



34 


the form 

= iO(h/Max. Asc. Rate) 

(3.8) 

7^2 = 0 . 

Similarly, for negative glide slope the penalty functions 
are 

T'l = 0 

(3.9) 

^2 ~ Des. Rate), 

The logic within the program has been designed to test h 
and determine which of the forms (3.8) or (3.9) is to be 
used. 

3.5,2) Altitude Constraints 


As previously mentioned, the aerodynamic model ignores 

the ground effect which is negligible above a relative 

altitude of 100 feet. Thus, the minimum altitude (h . ) 

min 

for which aircraft performance can be accurately simulated 
^^min- = ^ " \irnort^* penalty function which 


<Vn = h - ' 

enforces this condition is 

'f’3 = 

airport 


(3.. 10) 


The value is the mean sea level altitude of the 

airport and is specified in the program input data. 


*^he maximum altitude (h^^^^) depends upon the status 
of the flaps. If the flaps are not deployed, the limit 
is 35,000 feetj otherwise, it is 20, 000 feet. The penalty 
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function has the form 

Va = (3.11) 

max 

3.5,3) Load Factor 


The structural limits of the aircraft must be enforced 
to prevent trajectories which require forces and moments 
that damage the vehicle. Typically, this constraint is 
displayed as the load factor limits shown in Figures 3.7 


and 3.8. The load factor penalty for undeployed flaps is 




(3.12) 


Whenever the load factor is outside the' region -l<n <2.5, 

z 

the structural limit of the aircraft has been exceeded. 
When the flaps are deployed the penalty is 

■ = 10 |(n^ - 1.0)|. (3.13) 


3.5,4) Stall .and Buffet 


Figure 3,9 presents the stall speed characteristics of 
the .aircraft. These curves are approximated by the linear 
functions 

= 123.5 + 7.75 Kj flaps = 0 

= 99.25 + 6.75 Kj flaps = 1 (3.l4a) 

= 98,50 + 6.88 Kj flaps = 2 
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V 


= 97.00 + 6.63 K 
= 94.00 + 6.38 K 
= 92.00 + 6.25 K 
= 90.50 + 6.50 K 
= 86.25 + 6.63 K 
= 80.00 + 7.60 K 


flaps = 5 
flaps = 10 
flaps = 15 
flaps = 25 
flaps -- 30 
flaps - 45, 


(3.14b) 


■where is the stall speed in knots, W is the aircraft 
gross -weight in pouiids, and 


K = 10 W - 70,000) . 


(3.14c) 


The stall speed for intermediate flap settings is obtained 
by linear interpolation. These equations apply for -unde- 
ployed flight spoilers. If the spoilers are deployed, ' an 
additional 1 knot is added to the value calculated from 
(3.14). The stall speed was assumed to be independent of 
both the center of gravity and the thrust level. The penalty 
function was developed to maintain a ten percent safety mar- 
gin for V_. 

1.1 V 

fe = ) (3.15) 

This constraint is violated when the aircraft speed is 


lower than I.IV , 

s 

Figure 3.10 shows the low speed buffet angle of attack 
(ttg) as a function of the flap setting for undeployed 
spoilers. Linear interpolation between the discrete flap 
settings is used for intermediate values. The penalty func- 
tion for, low speed buffet is 

V^7 = 


(3.16) 
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If the flight spoilers are deployed, then three degrees are 
subtracted from the calculated value of (a„). This modified 
value is then used in (3.16). 

The high speed buffet constraint is shown in Figure 3.11 
to be a limitation on the lift coefficient. This maximum 
lift coefficient is approximated by a linear depen- 

dence on the Mach niimber of the form 

Clp = 1.342 - 0.7122 M. (3.17) 

This relationship applies for Mach numbers above 0.24. 

Below this value,. (3.17) does not apply and C^. is set 

lip 

equal to zero (yg=0). The penalty function for this con- 
straint is 

C. 

fQ " (3.18) 

Lp 

where is the current lift coefficient being tested. 

3.5.5) Angle of Attack 


An additional constraint on the angle of attack was ' 
formulated to insure the validity of the functions which 
represent the stability derivatives. A large proportion of 
the derivatives are functions of If this angle 

assumes values outside the range (-5°^a^(^p^25°) the deriva- 
tives calculated would be grossly in error. To prevent this 


possibility the penalty function 
) 

9 


= 10 — ) 


15) 


(3.19) 
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is used. 

3.6) Aircraft Control Surface Constraints 

The- aircraft control surfaces; i.e., elevator, rudder, 
etc., are subject to constraints which limit the. vehicle’s 
capability. These constraints are in the form of maximum 
displacements and displacement rates as described in this 
section. Table III presents the constraints for the vari- 
ous control surfaces under no load conditions. It was 
assumed that these rates also applied for loaded conditions. 


Table III 

Maximum Control Deflections and Actuation Rates 


Control 

Surface 

.J4aximum 

Displacement 

(deg) 

Maximum 

Rate 

(deg/sec) 

Elevator 

±21 

±56 

Stabilizer 

. 0-2.6 

. ±.56 - 

Manual Trim Wheel 

0-17.0 

±3.27 

Ailerons 

±20 

.‘±66 

(Panels 1,2, 3, 6, 7, 8) 
Spoilers 

0-40.0 

±60 

Rudder 

±24 

±56 
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The rate limit on the manual trim wheel cited in 

Table III is an estimated value. This estimate is based 

upon the assumption that the trim wheel displacement (s ) 

is proportional to the stabilizer deflection (6 ). 

s 

®p = ^ (3.20) 

The constant of proportionality, k, is assumed to be the 
ratio of the maximum displacements of the trim wheel to 
the stabilizer (k = 17/2.6 = 6.54). The trim wheel deflec- 
tion rate is then obtained as the time derivative of (3.20) 

« 

and noting that 6^ = 0.50 deg/sec in Table III. 

The flap constraints have been modeled in a more com- 
plex manner..' Table IV presents the flap displacement rate 
constraints. The operation times cited in the table are 
those times for which the flaps reached the particular posi- 
tion when started from an undeployed state at zero time. 

The relative operation times in Table IV are considered 
constant} i.e., the time for the flaps to go from a position 
of 5 to 10 is the same as from 10 to 5, which is 4.37 sec- 
onds. It is important to realize that the simulated flap 
position reflects the current stabilizer • deflection and not 
actual pilot guadrant setting . The linkage system 
connecting the pilot quadrant to the stabilizer is assumed 
to have a negligible delay time. 

The leading -edges of the flaps have different opera- 
tion times then those of Table IV. However, their contri- 
butions to the stability derivatives were not isolated in 
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Table IV 

Flap Operation Times 


Flap Position on Pilot Normal Operation 

Flap Quadrant Time, (sec) 


0 

o 

• 

o 

1 

5.20 

2 

10.12 

5 

21.67 

10 

26.04 

15 

28.65 

25 

29.80 

30 

32.00 

40 ' 

35.00 
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the data. Therefore, the assumption was made that all the 
flap aerodynamics were concentrated in the trailing edge. 
This allowed the leading edge operatioiiitimee to he ig- '• 
nored. . 

The program logic was designed to test each control 
vector element to determine if the magnitude and rate' con- 
straints are satisfied. If the constraints are violated 
at a particular instant, the control value is adjusted to 
the constraint boundary value. 

3.7) Aircraft Engine Constraints 

Two Pratt & Whitney JT8D-7 turbofan engines are mounted 
on the wings and are symmetrically positioned relative to 
the aircraft vertical plane of symmetry. The characteris- 
tics of these engines are presented in this section. 

3.7,1) Operational Limits 

The range of thrust per engine was assumed to be 1,540 
to 14,000 pounds. The simulation requires an accurate model 
of the engines because of the anticipated optimal trajec- 
tories. The optimization algorithm can be expected to 
generate flight paths which contain steep glide slopes. It 
is possible that the thrust requirements of an unconstrained 
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optimal trajectory ■will exceed the capability of the engines. 
Therefore, accurate models of the engine acceleration and 
deceleration characteristics are necessary. 

Figure 3.12 displays the engine deceleration character- 
istics. The curves beginning at thrust levels below sixty 
percent of the maximum rated thrust were extrapolated from 
the original curves . The curve originating at the one 
hundred percent level has been approximated by 

'i - 1.0 - 0.95 t + 0.4125 t^ - 0.0625 t^j0<t£i3 
. T* - < = 0.385 - 0.0925 t + 0.0075 t^;3<t^5 (3.2l) 

= 0.11;t>5, 

V 

t 

where T is the thrust per engine divided by 14,000 pounds. 
For the eighty percent level 

= 0.80 - 0.70 t + 0.295 t^ - 0.0442 t^;0£t^3 
t' < = 0.325 - 0.07 t + 0,005 t^;3£t*5 (3.22) 

' = 0.1l5t>5, 

for the sixty percent level 

= 0.60 - 0.4517 t + 0.1775 t^ - 0.0258 t^j0<t 3 
. T* / = 0.253 - 0.0405 t + 0.0015 t^;3<t*5- (3.23) 

- 0.11;'t>5, 

for the forty percent level 

■ I' = 0.40 - 0.1125 t - 0.015 t^ + 0.-0075 t^j0<t^3 
T* = 0.22 - 0.03t;3^t^4 
= 0.11;t>4, 


(3.24) 
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for the twenty percent level 

, j = 0.20 - 0.0633 t + 0.015 t^ - 0.001667 t^;0<ti:2.6 

■ V 

= 0.11jt>2.6, (3.25) 

and for thrust levels below twenty percent 

, V = T - 0.0633 t + 0.015 - 0.001667 t^;0<ts2.6 

T (. 

■I = 0.11jt>2.6 (3.26) 

where T is the thrust level coinciding with the initiation 

of the engine deceleration sequence. 

The program logic was constructed to identify the 

thrust level at the initiation of an engine deceleration 

interval. The thrust constraints for subsequent deceleration 

times are calculated using linear interpolation. The actual 

thrust level at each time is compared to the constraint 

magnitude. If the constraint is violated, the thrust is 

adjusted to that value of the constraint? otherwise it is 

left unchanged. In this manner the simulated thrust history 

is assured to have realistic deceleration characteristics. 

The engine acceleration curves are shown in Figure 3.13. 

The functional approximation of the curves originating at 

the ten percent level is 

■ / = 0.10 + t ; 0<t^l .,5 

T* <^ = 0.11 + y ^•^ ;1.5<t<2.85 (3.26) 

^ =5- 0.285 + (t - 2.85) 0.38;t>2,85, 
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for the twenty-five percent level 
= 0,25 + 0.01 t;0<t^l.0 
= 0.26 + 0.38(t-1.0) it>1.0_, 
for the fifty percent level 


- 0.50 + 


t j 0<t£0.65 


65.0 

= 0.51 + 0.38(t-0.65) ;t>0.65, 

for the seventy-five percent level 

1 


T 


(3.27) 


(3.28) 


= 0,75 + Q - tj0<t^0.45 

0.76 + 0,38(t-0.45) jt>0.45, 
fo,r an initial thrust "T above the seventy-five percent 
level 


(3.29) 


T=f+0.38tjt>0. (3.30) 

The thrust levels during an engine acceleration sequence 
are determined by the same method used for hhe deceleration 
case. In this manner a realistic simulation of the engine 
acceleration characteristics is achieved. 


3,^,2) Fuel Consumption Model 


The performance index contains a component which repre- 
sents the fuel consumed during a simulation. The perfor- 
mance index was formulated to permit arbitrary weighting.'of 
the various components (time, fuel, and noise). If the 
fuel consumption is a significant component an accurate 
model of the flowrate is necessary. Thrust specific fuel 
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consumption (TSFC) data, obtained from NASA Langley per- 
sonnel, are shown in Figures 3.14 through 3.16. [IQ],. -The. data 
displays the dependence of TSFC on Mach number, thrust, and 
altitude. These curves are strictly valid for the altitudes 
(•10,000; 5,000; 0). For altitudes above 10,000 feet the 
TSFC curves for the latter altitude were used. If the Mach 
number exceeds 0,6, the 0 . 6 Mach number data are used , 

For a given altitude and Mach number, the TSFC, . having' 
dimensions of Ibs/hr/lbf, was obtained as a function of the 
thrust using a least -squares polynomial criterion. For each 
discrete Mach number (0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6) and 
altitude (0, 5000, 10000 feet), a separate function is used. 
The approximation polynomials at sea level are 

= 1.244 - 0.3167 F + 0.0545 - 4.0 x lO”^ F^ 

+ 1.08 X 10“"^ f"^; M=0. 0,0.1 

= 1.390 - 0.3534 F + 0.0599 F^ - 4.385 x 10“^ F^ 

+ 1.18 X 10“"^ F^; M=0.2 

= 1.704 - 0.4973 F -f 0.0884 - 6.765 x lO"^ F^ 

+ 1.88 x.lO"^ F^; M=0.3 (3.31) 

= 1.896 - 0.5534 F + 0.0964 F^ - 7.25 x 10“^ F^ 

+ 1.99 X 10“^^ f"^; M=0.4 
= 2.653 - 1.137 F + 0..2784 F^ - 0.0334 F^ 

+ 1.933 X 10“^ f"^ - 4.306 x lO”^ F^j M=0,5 
= 4.58 - 2.76 F + 0.56 F^; M=0.6, F<2 . 5 
= 1.949 - 0.4641 F + 0.0729 F^ - 5.109 x 10“^ F^ 

+ 1.334 X 10"^ F'^; M=0.6, F>2.5, 



TSFC 



Net Thrust (lbs) 
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JT8D-7,7A Turbofan. Engine 
Standard Atmospheric Conditions 



Figure 3,14 

Thrust Specific Fuel Consumption at Sea Level 



Net Thrust (lbs) 
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JT8D-7,7A Turbofan Engine 



Figure 3.1c- 

Thrust Specific Fuel Consumption at 5,000 Feet 
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JT8D-7,7A Turbofan Engine 



Figure 3,15 


Thrust Specific Fuel Cons’rmption at 10,000 Feet 
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■where F is the thrust of one engine divided by 1,000 pounds. 
The approximation polynomials. for an altitude of 5,000 feet 


are 


TSFC 


L_= 1.133 - 0.29 If + 0.0527 - 4.05 x lO"^ F^ • 

+ 1.144 f"^j M=0. 0,0.1 

= 1.353 - 0.3769 F + 0.0685 F^ - 5.267 x lO"^ 

+ 1.48 X 10“^ F^j M=0.2 

= 1.453 - 0.3855 F + 0.0678 F^ - 5.11 x 10“^ F^ 
y + 1.42 X 10"^ F^j M=0.3 (3.32) 

■ = 1.568 - 0.3989 F + 0.06728 F^ - 4.89 x lO"^ 

+ 1.32 X 10 f"^; M=0.4 

= 2.174 - 0.8681 F + 0.2147 F^ - 0.0261 F^ 

+ 1.54 X 10"^ f"^ - 3.489 x 10“^ F^j M=0.5 
= 3.83 - 1.32 Fj M=0.6, F<2.0 

= 1.774 - 0.4202 F + 0.06649 F^ - 4.61 x 10“^ F^ 
■+ 1.19 X 10“^ F^; M=0.6, F>2.0, 


and an altitude of 10,000 feet are 


TSFC < 




1.062 - 0.2934 F + 0.06 F^ - 5.19 x 10“^ F^ 

+ 1.674 X 10"“^ F^j M=0. 0,0.1 

1.2329 ~ 0.3578 F + 0.0729 F^ - 6.356 x lO"^ F 
+ 2.06 X 10"^ F^j M=0.2 (3.33) 

1.2398 - 0.3038 F + 0.0569 F^ - 4.623 x lO*"^ F^ 
+ 1.428 X lO^^F^j M=0.3 

1.433 - 0.3868 F + 0.07376 F^ - 6.08 x 10“^ F^ 

+ 1.88 X lO”'^ f"^,* M=0.4 

1.575 - 0.4271 F + 0.0793 F^ - 6.38 x lO"^ F^ 

+ 1.92 X 10”^ F'^; M=0.5 
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TSFC = 1.7095 - 0.4778 F + 0.09 F^ - 7.376 x 10“^ F^ 

+ 2.24 X lO""^ F^; M=0.6. 

The TSFC is converted into fuel flowrate (ibs/sec) by the 
relationship 

m =: TSFC (2 engines) (Thrust/engine )/( 3, 600 sec/hr) (3.34) 
Interpolation is used for Mach numbers and altitude 
other than those for which (3.31) through (3.33) apply. 

As an example, assiame the value for TSFC is desired at a 
Mach number of 0.45, an altitude of 3,000 feet, and a thrust 
of 4,000 pounds. First the values of TSFC at sea level, 
for Mach numbers of 0.4 and 0.5 are calculated > to be 0.811 
and 0.873, respectively. The interpolated value at a Mach 
number of 0.45 and sea level altitude is 0.842. Similarly, 
the interpolated value at 5,000 feet is 0.784. Interpola- 
ting for the altitude of 3,000 feet and a Mach number of 0.45 
gives 0.807. In terms of the fuel flowrate, this is 1.79 
Ibs/sec. 

This model could be considered too complex for some 
9-PPlications of the optimization program. Specifically, 

.when the fuel component of the performance index is rela- 
tively small, the model is inefficient. However, the 
relative sizes of the components are determined by the pro- 
gram user, and the more accurate model presented above is 
necessary. 
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3.8) Model Accuracy 

The previous sections have detailed the development of 
a software model for a Boeing 737 aircraft. The foundation 
upon which the model has been built is the least-squares 
polynomial approximation of the stability derivatives. If 
this representation of the aerodynamic data does not realis- 
tically simulate the aircraft performance, then the gener- 
ated optimal trajectories are meaningless. Therefore, 
considerable effort was expended to verify the model's 
accuracy . 

Four cases have been chosen to evaluate the model accu- 
racy. Case #1 is the aircraft state at one instant along 
a -3 degree glide slope confined to the vertical plane. 

Case #2 is from a -6 degree glide slope which was also 
restricted to only longitudinal motion. Case #3 is from a 
more general trajectory confined to longitudinal motion. 
Finally, Case #4 is from a 20 degree banKed turn. The air- 
craft state and properties for each of the cases are summa- 
rized in Table. V. 

Table VI is a comparison of the lift coefficients of 
the model with a manual estimate obtained from the original 
data given in [sj. The results for each of the four cases 
demonstrate an excellent correspondence between the model 
and the actual values for the total lift coefficients. 
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Examining the component due to indicates a large error 

in the model’s representation of (<^Cj^/^a). However, the 
relatively small contribution which this component has 
toward_hhe total value limits the effect of the error to 

. * ’ t i 

an acceptable level. 

Tables VII and VIII compare the model and original 
data for the drag and lateral force coefficients, respec- 

I 

tively. The drag coefficient in the model has a large error 
for the sideslip component. The sideslip contributes approx' 
imately 0.03 percent of the total drag coefficient in case' 

#4, and a large error in representing this component is 
unimportant. The lateral force coefficient has all of its 
components accurately modeled. Therefore, the aerodynamic 
forces (lift, drag, and lateral) can be assumed to be accu- 
rately represented by the model. 

Tables IX, X, and XI show the comparison for the pitch, 
roTl, and yaw moment coefficients, . respectively. The worst 
aggregate error in the pitch model occurs for' case #4. The 
13.3 percent error is marginally acceptable. Currently this 
error will be tolerated, but if future circumstances warrant 
it, the model accuracy can be improved. The roll moment 
coefficient also has a large aggregate model error. However 
closer examination of Table X indicates that a near equi- 
librium condition of zero roll occurs. The absolute value 
of the 156 percent error can be estimated to 'cause 
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Comparison of the Drag Coefficients 
of the Model and the Original Data 
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Table IX 

Comparison of the Pitch Coefficients 
of the Model and the Original Data 
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Table XI 

Comparison of the Yaw Coefficients 
of the Model and the Original Data 
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a small angular acceleration. 

^^roll “ (rad/sGc^) (3.35) 

Assuming 

P = 2-. 38 X 10”^ slugs/ft? (sea level) 

V = 202 ft/sec 
S = 980 
b = 93 ft 

Cl = ( 1.229 - (-2,18)) X 10“"^ = 3.409 x 10”“^ 

= 375,000 slugs-ft^ 

then 

^^roll ~ ft-lb/ 375,000 slug-ft^ (3.36) 

or 

^ ^roll ~ X 10 rad/sec = 0,07 deg/sec^ (3.37) 

Thxs error in the acceleration is very small. Employing 
(3.35) -with the substitution of 

^zz “ 1»200,000 slugs-ft^ 

= (6.943 - 6.208) x lO"^ = 7.35 x lO"^ 

^xx *^ 1 * respectively, gives the yaw acceleration 

'^^yaw == ft-lb/ 1,200,000 slugs -ft^ (3.38) 

or 

^ ^yaw “ 2.7 X 10 rad/sec = 0.16 deg/sec^ (3.39) 
This error is also small. Therefore, the model accuracy 
for all three moment coefficients is acceptable. 
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3.9) Summary 

The aircraft model was developed utilizing least-squares 
polynomials approximating the stability derivatives for a 
Boeing 737, Fourier cosine series were also applied, but 
only for five derivatives. In addition to modeling the 
derivatives, the aircraft's state and control constraints 
were also included. The state constraints were enforced 
either by constants or piecewise polynomials. At each dis- 
crete time point the control vector was compared to the 
constraints and adjusted to insure compliance, if necessary. 
Both magnitude and rate constraints on the controls were 
considered. 

A ntanber of checks were conducted on the model to 
verify its correspondence to the actual aircraft response 
characteristics. In all of these tests, the calculated aero- 
dynamic forces and moments in the model closely paralleled 
those of the actual aircraft. 



IV. 


PROGRAMMING THE STEEPEST ■ 


DESCENT OPTIMIZATION PROCEDURE 


4 . 1 Introduction 

The annual report [2] submitted a year ago summarized 
the mathemetics involved in the steepest descent procedure. 
The equations now have all been programmed and apparently 
debugged. 

4 . 2 Problem Areas 

During the first attempts to run the program, several 
problems arose. One was related to the passenger comfort 
model. A measure of passenger discomfort had been normalized 
to unity and then raised to the power of one hundred. If 
there was passenger discomfort, then the indicator would 
exceed unity and, when raised to the power of one hundred, 
would be extremely large. Likewise, when the indicator was 
less than unity, the quantity raised to the power of one hun- 
dred would be extremely small. This type behavior was 
thought to be desirable for use as a penalty function to be 
included with the performance index; however, when there was 
discomfort, the indicator raised to the power of one hundred 
became so large that it exceeded the handling capacity of the 
computer. It was then decided to raise the indicator to the 
power of ten if the indicator exceeded unity and to the power 
of one hundred if the indicator was less than unity. This 
seemed to work and still accomplish the objective. 
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Another area in which we have experienced some diffi- 
culty has to do with the stopping condition. Ideally, the 
optimization procedure should be allowed to pick that termi- 
nal time which is optimum for the situation. To allow this 
freedom a stopping condition must be chosen. As was des- 
cribed in the final report [2 for last year's work, the 
stopping condition was formulated as the time rate change of 
distance between the aircraft and the airport. Whenever this 
distance stopped decreasing and began increasing (i.e. when 
its time rate of change went through zero) , the trajectory 
would terminate. Theoretically, successive iterations would 
cause this event to occur closer and closer to the desired 
boundary conditions. In order to prevent over-restricting 
the maneuver, it was decided to begin testing for the stop- 
ping condition only after the aircraft had gotten within five 

miles of the airport. 

The net result was that the aircraft seemed to seek out 
areas of low population density and circle over these areas. 
This would continue until the upper limit on computer time 
was reached at which time the trajectory was terminated. 
Attempts to resolve this problem by additional emphasis on- 
boundary condition errors and time helped somewhat but did 
not cure the problem. It was then decided to switch over to 
a fixed-time approach. The time required for a straight-in 
constant-glide-slope trajectory was calculated. This time 
was increased by twenty percent to allow time for curves in 
the trajectory. The problem was then run as one of fixed time. 



4.3 Program Size and Cost 


As one would guess, the program for performing the opti- 
mization is quite complex and lengthy. The total storage 
requirements vary, depending on the length of the trajectory 
to be optimized. For a 500-secpnd trajectory (20 miles) the 

3 

storage requirements are 41 x 10 words. The total computa- 
tion time per iteration for a 500-second trajectory is 750 
seconds. Note that this includes not only the aircraft equa- 
tions of motion but the evaluation of the noise effect every 
five seconds and also the backward integration of three sets 
of adjoint equations. The total computing cost per iteration 
for a 500-second trajectory is $40. We presently estimate 
that twenty -iterations may be required to converge to the 
optimal solution. This comes to $800 per optimal trajectory. 
The entry point into the near terminal area is variable, and 
there are four runways which can be used for landing. Thus, 
optimal trajectories for several sets of boundary conditions 
need to be calculated. The next section discusses the 
results obtained so far. 



V. RESULTS 


5.1 Introduction 

The utility of a performance index for comparing various 
trajectories is apparent and has been discussed in some 
detail in the report [2] of our work in 1974-75. There, we 
reported comparisons of three-degree glide-slope trajectories 
with six-degree glide-slope trajectories. 

Since that time, we have been able to implement the 
optimization procedure and allow it to search for the optimal 
trajectory. We began with a one-hundred second trajectory 
(approximately four miles) . The initial trajectory was 
straight in at a glide slope of three degrees. Table XII 
lists the results of ten consecutive iterations. 

The quantity J is the performance measure. It is given 
by the equation 

J = . 7 X Time (sec) + . 05 x Fuel (lbs) + .0001 x Noise (People-sec) 

The next column lists .25 dij;, where diff is the error in 
boundary conditions. This consists of 

di> = .25x(X - X^)^ + .25x(Y - Y^) ^ + .25x(Z - Z^) ^ 

+ 625x(/- y'g)^ + 5000x(i{i - ^^)^ 

where X, Y, and Z represent the coordinates of the aircraft 
at final time, and X^, Y^, and Z^ represent the desired 
values for these coordinates. The angle Y is the flight path 
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ITERATION 

J 

0.25 

X 

di]^ 

A 

■^PRED 

Xp(ft) 

-6976 

Yp(ft) 

5654 


YpCdeg) 

-2.75 

li)p(rad) 

-0.60 

1 

103.0 

..2947 

X 

10^ 

+ 

7.62 

-5,248 

4,476 

993 

-2.41 

- .60 

2 

107.6 

.7795 

X 

10^ 

+ 

0.25 

' -16,539 

- 91.6 

815 

-5.41 

-2.50 

3 

108.2 

,3979 

X 

10^ 

- 

0.52 

-12,222 

- 343 

790 

-5.52 

-1.67 

4 

106.9 

.2324 

X 

10^ 

- 

1.37 

' -10,013 

380 

739 

-5.45 

-1.29 

5 

103.7 

.1341 

X 

10^ 

- 

1.98 

- 8,812 

1,414 

695 

-5.45 

-1.38 

6 

99.7 

.7498 

X 

10^ 

- 

2.21 

- 8,244 

2,446 

651 

-5.48 

-1.61 

7 

96.6 

.4307 

X 

10^ 

- 

1.30 

- 8,320 

3,418 

592 

-5.32 

-1.91 

8 

94.8 

.2698 

X 

10^ 


0.85 

- 8,509 

4,278 

502 

-5.10 

-2.16 

9 

93.3 

.1620 

X 

10^ 

- 

0.83 

- 8,388 

4,923 

402 

-4.61 

-2.26 

10 

92.0 

. 14089 

X 

10^ 


— 

- 8,428 

5,896 

289 

-4.14 

-2.41 


Table XII 

Results of Ten Iterations on A Four-Mile Trajectory 
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angle, and ti» is the heading angle. 

5 . 2 Pattern of Convergence 

One important feature to note is that J the performance 
measure decreased monitomically after the third iteration. 
Also, the boundary condition error, dij) decreased monitorai- 
caliy after the second iteration. The predicted change in 
performance, always had the correct sign and, in most 

cases, was quite close in magnitude. This closeness is a 
measure of the degree with which the linearity assumptions 
are being met. If these assumptions are grossly violated, 
there is no guarantee of convergence toward an optimum. To 
ensure linearity one must keep the step size from one itera- 
tion to the next sufficiently small. Our program automati- 
cally reduces the allowable step size to one half its previ- 
ous value if the performance measure ever increases. ■ So far, 
■this strategy seems to be working satisfactorily. 

'Examining once more the column on performance index, J 
.it appears that reducing the value from 103 to 92 is not all 
that significant. However, recall that we are now working 
with a fixed time problem which means that a portion of J , 

.7 X 104 seconds, is fixed. The variable portion of the 
performance measure actually decreased from 30 on the first 
iteration to 19 on the tenth iteration. The noise component 
decreased from 21 to 11. Thus, a significant improvement 
has been achieved. 

Judging from dt(), it is apparent that the convergence is. 
not complete, i.e. the boundary conditions have not been 
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completely satisfied. Figure 5.1 illustrates this also. 
However, it is believed that additional iterations would 
bring these boundary condition errors closer to zero. On 
Fig. 5.1 it is seen that the trajectory for the tenth itera- 
tion is one which swings across the areas of low population 
density and cuts between the two areas of high population 
density. This is as one would hope the procedure would work. 
Hopefully, additional iterations would change the curvature 
at. the very end of the trajectory and make the final heading 
of the aircraft toward the runway. 

The four-mile trajectory is actually not a realistic 
problem, although it does seem to prove out the workability 
of our program. Such a short trajectory does not allow room 
for very much maneuvering or time to recover from the maneu- 
•ver and align with the runway. By going to the twenty-mile 
trajectory, more freedom will be given to the procedure? and 
.the areas of low population density may be sought out while 
still affording time for the aircraft to straighten out for 
the final approach to the airport. A shortage of computing 
funds precluded the optimization of any twenty-mile trajec- 
tories this past year. 



Results of Ten 
on a Four-Mile 






VI. CONCLUSIONS AND FUTURE PLANS 


At this point all our models have been developed and 
programmed, and our procedure is working. There is still a 
certain amount of art as well as science involved in the suc- 
cessful use of the procedure. Convergence is very much a 
function of the initial guess and the step size. Further- 
more, a step size which works for awhile may cause oscilla- 
tion later and have to be reduced. In some cases automatic 
reduction of step size has been successful; however, in gen- 
eral a certain amount of user interaction seems to be a 
necessity. 

The results we have obtained for the four— mile trajec- 
tory are very encouraging. Things seem to be behaving as one 
would expect. We are most anxious to begin applying our pro- 
gram to the twenty-mile trajectories. There are several sets 
of boundary conditions for which we desire optimal trajector- 
ies. If possible, we would also like to investigate the sen- 
sitivity of the optimal trajectories to changes in the values 
of the weighting parameters in the performance index. The 
inclusion of wind in our simulation is a possibility. One 
other item which may be worth considering is the frequency of 
landings as opposed to treating each landing separately. 
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